F105<- read.table("C:/Users/lixia/OneDrive/Documents/mémoire/Data/new/Female105.txt", header=TRUE, sep="\t", dec=",")

M105<- read.table("C:/Users/lixia/OneDrive/Documents/mémoire/Data/new/Male105.txt", header=TRUE, sep="\t", dec=",")

attach(F105)
attach(M105)


## remove the NA
nona_BEL_F = BEL_F[!is.na(BEL_F)]
nona_DEU_F = DEU_F[!is.na(DEU_F)]
nona_CHE_F = CHE_F[!is.na(CHE_F)]
nona_CAN_F = CAN_F[!is.na(CAN_F)]
nona_USA_F = USA_F[!is.na(USA_F)]


nona_BEL_M = BEL_M[!is.na(BEL_M)]
nona_DEU_M = DEU_M[!is.na(DEU_M)]
nona_CHE_M = CHE_M[!is.na(CHE_M)]
nona_CAN_M = CAN_M[!is.na(CAN_M)]
nona_USA_M = USA_M[!is.na(USA_M)]



############# K-sample AD Tests to verify if the data of different countries is from the same unknown distribution################

library("kSamples")

###Female###
ad.test(nona_DEU_F,nona_BEL_F,nona_CHE_F, FRA_F, method = "asymptotic", dist = FALSE)
##p-value: 0.02739 !!! 



#why##
ks.test(nona_DEU_F, FRA_F, 
        alternative = "less",
        exact = NULL)
ad.test(nona_DEU_F,nona_BEL_F,nona_CHE_F, method = "asymptotic", dist = TRUE)
##p-value: 0.6911
ad.test(nona_BEL_F,nona_CHE_F, FRA_F, method = "asymptotic", dist = FALSE)
##p-value: 0.6042 

summary(FRA_F)
summary(nona_DEU_F)


###Male###
ad.test(nona_DEU_M,nona_BEL_M,nona_CHE_M,FRA_M, method = "asymptotic", dist = FALSE)
##P-value: 0.4664




##Histogram

hist(nona_BEL_F, breaks=40, 
     freq=FALSE, xlim=c(105, 115),
     main="BE_Female",
     xlab="x")
hist(nona_DEU_F, breaks=40, 
     freq=FALSE, xlim=c(105, 115),
     main="DEU_Female",
     xlab="x")
hist(FRA_F, breaks=40, 
     freq=FALSE, xlim=c(105, 117),
     main="FR_Female",
     xlab="x")

hist(nona_CHE_F,breaks=40, 
     freq=FALSE, xlim=c(105, 115),
     main="CHE_Female",
     xlab="x")

##hist(nona_CAN_F,breaks=40, 
     freq=FALSE, xlim=c(105, 116),
     main="CAN_Female",
     xlab="x")



hist(nona_BEL_M, breaks=10, 
     freq=FALSE, xlim=c(105, 115),
     main="BE_Male",
     xlab="x")
hist(nona_DEU_M, breaks=20, 
     freq=FALSE, xlim=c(105, 115),
     main="DEU_Male",
     xlab="x")
hist(FRA_M, breaks=20, 
     freq=FALSE, xlim=c(105, 117),
     main="FR_Male",
     xlab="x")

hist(nona_CHE_M,breaks=10, 
     freq=FALSE, xlim=c(105, 115),
     main="CHE_Male",
     xlab="x")

##hist(nona_CAN_M,breaks=10, 
     freq=FALSE, xlim=c(105, 116),
     main="CAN_Male",
     xlab="x")





##group BEL FRA DEU CHE under EU, 
EU<- read.table("C:/Users/lixia/OneDrive/Documents/mémoire/Data/new/FM105_EU.txt", header=TRUE, sep="\t", dec=",")
attach(EU)
names(EU)

boxplot(EU)
nona_Dage_M = Dage_M[!is.na(Dage_M)]


####Trend test example######
##library("trend")



######Threshold selection#########

####Method 1#####

##mean excess plot
library("POT")
mrlplot(Dage_F, main="Female_BEL_FRA_CHE_DEU")
mrlplot(Dage_M, main="Male_BEL_FRA_CHE_DEU")


##QQ plot##
library("car")
qqPlot(Dage_F, distribution="exp",main="Female_BEL_FRA_CHE_DEU")
qqPlot(Dage_M, distribution="exp",main="Male_BEL_FRA_CHE_DEU")

qqPlot(nona_CAN_F, distribution="exp", main="Female_CAN")
qqPlot(nona_CAN_M, distribution="exp", main="Male_CAN")



##Shape Plot
shape(Dage_F, start = 11, models =100, end = 3000, reverse=FALSE)
shape(nona_Dage_M , start = 11, models =80, reverse=FALSE)



####Method 2 CVM & AD raw P-value goodness-of-fit test######
library("eva")
(threshold <- seq(105,109,by=0.01))

##CVM
gpdSeqTests(Dage_F, thresholds=threshold, method="cvm") #105.79
gpdSeqTests(nona_Dage_M, thresholds=threshold, method="cvm")#105
##AD
gpdSeqTests(Dage_F, thresholds=threshold, method="ad")#105.89
gpdSeqTests(nona_Dage_M, thresholds=threshold, method="ad")#105.01




####Method 3 Automated selection: CVM&AD + Forward Stop rule ######
install.packages("selectiveInference")
library("selectiveInference")

#########################################################

####Female Group####

##CVM##
threshold_F_105 <- seq(105,113, by=0.01)
out_cvm_F_105 <- gpdSeqTests(Dage_F, thresholds=threshold_F_105, method="cvm")
out_cvm_F_105
forwardStop(out_cvm_F_105$p.values,alpha=0.1) #92, 105.91 2373, 0.7446

threshold_F_110 <- seq(110,113, by=0.01)
out_cvm_F_110 <- gpdSeqTests(Dage_F, thresholds=threshold_F_110, method="cvm")
out_cvm_F_110
forwardStop(out_cvm_F_110$p.values,alpha=0.1) #11, 110.10 109   0.311896



##AD##
out_ad_F_105 <- gpdSeqTests(Dage_F, thresholds=threshold_F_105, method="ad")
out_ad_F_105
forwardStop(out_ad_F_105$p.values,alpha=0.1)  #187   106.86  1238   0.3856

out_ad_F_110 <- gpdSeqTests(Dage_F, thresholds=threshold_F_110, method="ad")
out_ad_F_110
forwardStop(out_ad_F_110$p.values,alpha=0.1) #13 110.12   105   0.44027


##French subsets##
out_ad_FRA_F_110<- gpdSeqTests(FRA_F, thresholds=threshold_F_110, method="ad")
out_ad_FRA_F_110
forwardStop(out_ad_FRA_F_110$p.values,alpha=0.1) #15  110.14   85   0.3288


#####Male Group#####
threshold_M_105<-seq(105,110,by=0.01)
out_cvm_M_105 <- gpdSeqTests(nona_Dage_M, thresholds=threshold_M_105, method="cvm")
out_cvm_M_105
forwardStop(out_cvm_M_105$p.values,alpha=0.1)  #0  105   430   0.6245

threshold_M_107<-seq(107,110,by=0.01)
out_cvm_M_107 <- gpdSeqTests(nona_Dage_M, thresholds=threshold_M_107, method="cvm")
out_cvm_M_107
forwardStop(out_cvm_M_107$p.values,alpha=0.1)  #0   107  89   0.8458

threshold_M_108<-seq(108,110,by=0.01)
out_cvm_M_108 <- gpdSeqTests(nona_Dage_M, thresholds=threshold_M_108, method="cvm")
out_cvm_M_108
forwardStop(out_cvm_M_108$p.values,alpha=0.1)  #2   108.01   46  0.0974

threshold_M_109<-seq(109,110,by=0.01)
out_cvm_M_109 <- gpdSeqTests(nona_Dage_M, thresholds=threshold_M_109, method="cvm")
out_cvm_M_109
forwardStop(out_cvm_M_109$p.values,alpha=0.1)   #0   109     20   0.4027



###Estimate of parameters of GPD####

library("POT")
?fitgpd
mle_F<-fitgpd(Dage_F, 106.86, "mle")$param
mom_F<-fitgpd(Dage_F,106.86, "moments")$param
print(rbind(mom_F, mle_F))

mle_M <- fitgpd(nona_Dage_M, 105, "mle")$param
mom_M <- fitgpd(nona_Dage_M, 105, "moments")$param
print(rbind(mom_M, mle_M))


##Condidence Intervals
mle_F<-fitgpd(Dage_F, 106.86, "mle")
gpd.fiscale(mle_F, conf=0.95)  
gpd.pfshape(mle_F, conf=0.95)


############Model checking############"
f_F_MLE <- gpdFit(Dage_F, nextremes= 1238, method="mle")
plot(f_F_MLE)

f_M_MLE <- gpdFit(nona_Dage_M, nextremes= 431, method="mle")
plot(f_M_MLE)2


#OR#
Exceedance_F<- read.table("C:/Users/lixia/OneDrive/Documents/mémoire/Data/Exceedance_F.txt", header=TRUE, sep="\t", dec=",")
gpd_F <- fitdist(Exceedance_F$Exceedance_F, "gpd", start=list(shape=-0.04, scale=1.4293))
summary(gpd_F )
plot(gpd_F )
gofstat(gpd_F)

Exceedance_M<- read.table("C:/Users/lixia/OneDrive/Documents/mémoire/Data/Exceedance_M.txt", header=TRUE, sep="\t", dec=",")
gpd_M <- fitdist(Exceedance_M$Exceedance_M, "gpd", start=list(shape=-0.06, scale=1.3963))
summary(gpd_M )
plot(gpd_M )
gofstat(gpd_M)



################Subsets analysis############


####mean excess plot####
library("POT")
mrlplot(nona_BEL_F, main="Female_BEL")
mrlplot(nona_CHE_F, main="Female_CHE")
mrlplot(nona_DEU_F, main="Female_DEU")
mrlplot(FRA_F, main="Female_FRA")
par(mfrow=c(2,2))

shape(nona_BEL_F, models =50, reverse=FALSE, main="BE_F")
shape(nona_CHE_F, models =50, reverse=FALSE,main="CHE_F")
shape(nona_DEU_F, models =50, reverse=FALSE,main="DEU_F")
shape(FRA_F, models =100, reverse=FALSE,main="FRA_F")
par(mfrow=c(2,2))
?shape



##threshold selection

##Female subsets###
(threshold <- seq(105,109,by=0.01))

out_ad_DEU_F<- gpdSeqTests(nona_DEU_F, thresholds=threshold, method="ad")
out_ad_DEU_F
forwardStop(out_ad_DEU_F$p.values,alpha=0.1)

out_ad_BE_F<- gpdSeqTests(nona_BEL_F, thresholds=threshold, method="ad")
out_ad_BE_F
forwardStop(out_ad_BE_F$p.values,alpha=0.1)

out_ad_CHE_F<- gpdSeqTests(nona_CHE_F, thresholds=threshold, method="ad")
out_ad_CHE_F
forwardStop(out_ad_CHE_F$p.values,alpha=0.1)

out_ad_FRA_F<- gpdSeqTests(FRA_F, thresholds=threshold, method="ad")
out_ad_FRA_F
forwardStop(out_ad_FRA_F$p.values,alpha=0.1)


##### parameters estimation for GPD FIT####
(mle_BE<-fitgpd(nona_BEL_F, 105, "mle"))
(mle_FR<-fitgpd(FRA_F, 106.21, "mle"))
(mle_DEU<-fitgpd(nona_DEU_F, 105.15, "mle"))
(mle_CHE<-fitgpd(nona_CHE_F, 105.04, "mle"))


(BE_F_MLE <- gpdFit(nona_BEL_F, nextremes= 269, method="mle"))



#################MISC



library("KSgeneral")
ks.test(nona_BEL_F, FRA_F) ## p-value: 0.4819  D = 0.053547,
ks.test(nona_BEL_F, nona_DEU_F) ## p-value: 0.4448  D = 0.060548,
ks.test(nona_BEL_F, nona_CHE_F) ## p-value: 0.8345  D = 0.061462,
ks.test(FRA_F, nona_CHE_F)## p-value: 0.3425 D = 0.075123
ks.test(FRA_F, nona_DEU_F) ## p-value: 0.04479 D = 0.054292 !!!


ks.test(c(nona_DEU_F,nona_BEL_F,nona_CHE_F),FRA_F)## p-value: 0.05509
ks.test(c(nona_BEL_F,nona_CHE_F,FRA_F),nona_DEU_F) ## p-value = 0.05293

## Kolmogorov-Smirnov Tests to verify if the data among the birth years is homogenity
## Group by 6 years (group1: 1881-1886; group2: 1887-1892; group3: 1893-1898)

h_F <- read.table("C:/Users/lixia/OneDrive/Documents/mémoire/Data/new/homoge_year_F.txt", header=TRUE, sep="\t", dec=",")
attach(h_F)
names(h_F)
nona_Group.1_F = Group.1_F[!is.na(Group.1_F)]
nona_Group.2_F = Group.2_F[!is.na(Group.2_F)]
nona_Group.3_F = Group.3_F[!is.na(Group.3_F)]
ks.test(nona_Group.1_F, nona_Group.2_F) ## p-value: 0.04202
ks.test(nona_Group.1_F, nona_Group.3_F) ## p-value: 0.3083
ks.test(nona_Group.2_F, nona_Group.3_F) ## p-value: 0.3019

## group by 3 years
h1_F <- read.table("C:/Users/lixia/OneDrive/Documents/mémoire/Data/new/homoge_year_F1.txt", header=TRUE, sep="\t", dec=",")
attach(h1_F)
names(h1_F)
ks.test(G1_F, G2_F) ## p-value: 0.2389
ks.test(G1_F, G3_F) ## p-value: 0.01144
ks.test(G1_F, G4_F) ## p-value: 0.04887
ks.test(G1_F, G5_F) ## p-value: 0.09864
ks.test(G1_F, G6_F) ## p-value: 0.06407
ks.test(G2_F, G3_F)## p-value: 0.5041




